library('tidyverse')
library('tidyquant')
library('tseries')
library('gghighlight')
load('api_key.rdata')

Time series analysis:

What is a stationary time series?

“A time series is stationary if a shift in time doesn’t cause a change in the shape of the distribution. Basic properties of the distribution like the mean, variance, and covariance are constant over time.”

Stationary vs Seasonal Time series - http://www.statisticshowto.com/stationarity/

Stationary vs Seasonal Time series - http://www.statisticshowto.com/stationarity/

Different types of stationarity:

Why do we care about stationarity?

Most statistical models that have been developed assume stationarity.

set.seed(100)
# monthly data, normally distributed
# 100 stoltzcoin average trade price per month
# Standard deviation of 3
dates = seq(from = lubridate::ymd('2011-07-01'), 
            to = lubridate::ymd('2018-06-30'), 
            by = 'months')
n = length(dates)
mu = 100
sd = 3
time_series = tibble(date = dates, stoltzcoin = rnorm(n = n, mean = mu, sd = sd), 
                     type = 'stationary', trend = 'none')
ggplot(time_series, aes(x = date, y = stoltzcoin)) + 
  geom_line() + geom_smooth(method = lm, 
                            formula = y ~ x, se = FALSE) + 
  geom_hline(aes(yintercept = mean(stoltzcoin)), color="red")

Perform regression and look at the residuals

ts_fit = lm(stoltzcoin ~ date, data = time_series)
ts_fit

Call:
lm(formula = stoltzcoin ~ date, data = time_series)

Coefficients:
(Intercept)         date  
  1.082e+02   -5.002e-04  
res = tibble(residuals = resid(ts_fit))
ggplot(res, aes(x = 1:nrow(res), y = residuals)) + geom_point() + geom_rug(sides = 'l')

ggplot(res, aes(x = residuals)) + geom_histogram(aes(y = ..density..), bins = 40) + geom_density(col = 'red', size = 1.5)

Always test for stationarity even if it looks stationary

Augmented Dickey-Fuller is a typical test - tseries::adf.test()
  • Tests null hypothesis that the autoregressive model has a unit root
  • Stationarity means rejecting the null hypothesis
  • We will look for p-value < 0.05
# adf.test only accepts univariate time series data, zoo object will be used
adf.test(zoo::zoo(time_series$stoltzcoin, time_series$date))
p-value smaller than printed p-value

    Augmented Dickey-Fuller Test

data:  zoo::zoo(time_series$stoltzcoin, time_series$date)
Dickey-Fuller = -4.1967, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary

Add a linear trend for the example

Stoltzcoin has been going up by 2 each month

Known as “trend stationary”

going_up_by = 2
trend = seq(from = 0, 
            to = (n-1) * going_up_by,
            by = going_up_by)
ts_lin_trend = tibble(date = dates, stoltzcoin = time_series$stoltzcoin + trend, 
                     type = 'non-stationary', trend = 'linear')
ggplot(ts_lin_trend, aes(x = date, y = stoltzcoin)) + 
  geom_line() + geom_smooth(method = lm, 
                            formula = y ~ x, se = FALSE) + 
  geom_hline(aes(yintercept = mean(stoltzcoin)), color="red")

Perform regression and look at the residuals

ts_fit = lm(stoltzcoin ~ date, data = ts_lin_trend)
ts_fit

Call:
lm(formula = stoltzcoin ~ date, data = ts_lin_trend)

Coefficients:
(Intercept)         date  
 -887.75051      0.06521  
res = tibble(residuals = resid(ts_fit))
ggplot(res, aes(x = 1:nrow(res), y = residuals)) + geom_point() + geom_rug(sides = 'l')

ggplot(res, aes(x = residuals)) + geom_histogram(aes(y = ..density..), bins = 40) + geom_density(col = 'red', size = 1.5)

adf.test(zoo::zoo(ts_lin_trend$stoltzcoin, ts_lin_trend$date))
p-value smaller than printed p-value

    Augmented Dickey-Fuller Test

data:  zoo::zoo(ts_lin_trend$stoltzcoin, ts_lin_trend$date)
Dickey-Fuller = -4.1967, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary

Add a seasonal component to the linear trend

Stoltzcoin fluctuates throughout the year

going_up_by = 2
trend = seq(from = 0, 
            to = (n-1) * going_up_by,
            by = going_up_by)
seasonal = sin(1:n) * 10
ts_lin_trend_seasonal = tibble(date = dates, stoltzcoin = time_series$stoltzcoin + trend + seasonal, 
                     type = 'non-stationary', trend = 'linear')
ggplot(ts_lin_trend_seasonal, aes(x = date, y = stoltzcoin)) + 
  geom_line() + geom_smooth(method = lm, 
                            formula = y ~ x, se = FALSE) + 
  geom_hline(aes(yintercept = mean(stoltzcoin)), color="red")

Perform regression and look at the residuals

ts_fit = lm(stoltzcoin ~ date, data = ts_lin_trend_seasonal)
ts_fit

Call:
lm(formula = stoltzcoin ~ date, data = ts_lin_trend_seasonal)

Coefficients:
(Intercept)         date  
 -887.83081      0.06523  
res = tibble(residuals = resid(ts_fit))
ggplot(res, aes(x = 1:nrow(res), y = residuals)) + geom_point() + geom_rug(sides = 'l')

ggplot(res, aes(x = residuals)) + geom_histogram(aes(y = ..density..), bins = 40) + geom_density(col = 'red', size = 1.5)

adf.test(zoo::zoo(ts_lin_trend_seasonal$stoltzcoin, ts_lin_trend_seasonal$date))
p-value smaller than printed p-value

    Augmented Dickey-Fuller Test

data:  zoo::zoo(ts_lin_trend_seasonal$stoltzcoin, ts_lin_trend_seasonal$date)
Dickey-Fuller = -6.4213, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary

Add an exponential trend to the seasonal component

Stoltzcoin has increased in volatility over the years

  • Stoltzcoin now has:
    • A positive and growing linear trend
    • Seasonality
going_up_by = 2
trend = seq(from = 0, 
            to = (n-1) * going_up_by,
            by = going_up_by)
seasonal = sin(1:n) * 1:n
ts_lin_trend_seasonal_exp = tibble(date = dates, stoltzcoin = time_series$stoltzcoin + trend + seasonal, 
                     type = 'non-stationary', trend = 'linear')
ggplot(ts_lin_trend_seasonal_exp, aes(x = date, y = stoltzcoin)) + 
  geom_line() + geom_smooth(method = lm, 
                            formula = y ~ x, se = FALSE) + 
  geom_hline(aes(yintercept = mean(stoltzcoin)), color="red")

Perform regression and look at the residuals

ts_fit = lm(stoltzcoin ~ date, data = ts_lin_trend_seasonal_exp)
ts_fit

Call:
lm(formula = stoltzcoin ~ date, data = ts_lin_trend_seasonal_exp)

Coefficients:
(Intercept)         date  
 -925.35762      0.06756  
res = tibble(residuals = resid(ts_fit))
ggplot(res, aes(x = 1:nrow(res), y = residuals)) + geom_point() + geom_rug(sides = 'l')

ggplot(res, aes(x = residuals)) + geom_histogram(aes(y = ..density..), bins = 40) + geom_density(col = 'red', size = 1.5)

adf.test(zoo::zoo(ts_lin_trend_seasonal_exp$stoltzcoin, ts_lin_trend_seasonal_exp$date))
p-value smaller than printed p-value

    Augmented Dickey-Fuller Test

data:  zoo::zoo(ts_lin_trend_seasonal_exp$stoltzcoin, ts_lin_trend_seasonal_exp$date)
Dickey-Fuller = -6.2925, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary

R’s decompose function does a lot for you

d1y =  as.numeric(lubridate::year(min(ts_lin_trend_seasonal_exp$date)))
d1m = as.numeric(lubridate::month(min(ts_lin_trend_seasonal_exp$date)))
d2y =  as.numeric(lubridate::year(max(ts_lin_trend_seasonal_exp$date)))
d2m = as.numeric(lubridate::month(max(ts_lin_trend_seasonal_exp$date)))
x = ts(ts_lin_trend_seasonal_exp$stoltzcoin, start=c(d1y, d1m), end=c(d2y, d2m), frequency=12)
plot(decompose(x, type = 'additive'))

plot(decompose(x, type = 'multiplicative'))

Example from decompose function

random portion should just look like white noise (otherwise you’ve likely missed something)

require(graphics)
m <- decompose(co2)
plot(m)

dow = tq_index("DOW") %>%
  tq_get(get = "stock.prices")
head(dow, 10)
p = ggplot(dow, aes(x = date, y = close, col = symbol)) + 
  geom_line() + 
  gghighlight(max(close) > 250) + 
  theme_minimal() + 
  facet_wrap(~ symbol)
print(p)

Is BA a stationary time series?

Can you simply zoom in and find a portion that is stationary?

ba = dow %>% 
  filter(symbol == 'BA') %>% 
  select(date, close)
plotly::ggplotly(ggplot(ba, aes(x = date, y = close)) + geom_line() + geom_smooth(method = lm, se = FALSE))
ts_fit

Call:
lm(formula = close ~ date, data = ba)

Coefficients:
(Intercept)         date  
 -734.75085      0.05385  

Are the daily returns stationary?

ret = ba %>%
  tq_transmute(close, mutate_fun = dailyReturn)
ggplot(ret, aes(x = date, y = daily.returns)) + geom_line() + geom_smooth(method = lm, se = FALSE)

adf.test(zoo::zoo(ret$daily.returns, ret$date))
p-value smaller than printed p-value

    Augmented Dickey-Fuller Test

data:  zoo::zoo(ret$daily.returns, ret$date)
Dickey-Fuller = -14.308, Lag order = 13, p-value = 0.01
alternative hypothesis: stationary
ts_fit = lm(daily.returns ~ date, data = ret)
ts_fit

Call:
lm(formula = daily.returns ~ date, data = ret)

Coefficients:
(Intercept)         date  
 -8.376e-03    5.736e-07  
res = tibble(residuals = resid(ts_fit))
ggplot(res, aes(x = 1:nrow(res), y = residuals)) + geom_point() + geom_rug(sides = 'l')

ggplot(res, aes(x = residuals)) + geom_histogram(aes(y = ..density..), bins = 40) + geom_density(col = 'red', size = 1.5)

LS0tCnRpdGxlOiAiVGltZSBTZXJpZXMgTW9kZWxpbmcsIEFuYWx5c2lzIGFuZCBEZWNvbXBvc2l0aW9uIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpgYGB7ciwgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRX0KbGlicmFyeSgndGlkeXZlcnNlJykKbGlicmFyeSgndGlkeXF1YW50JykKbGlicmFyeSgndHNlcmllcycpCmxpYnJhcnkoJ2dnaGlnaGxpZ2h0JykKbG9hZCgnYXBpX2tleS5yZGF0YScpCmBgYAoKCiMjIyMgVGltZSBzZXJpZXMgYW5hbHlzaXM6ICAKCiAgKiBTdGF0aW9uYXJ5IHZzLiBOb24tU3RhdGlvbmFyeQogICogQWRkaXRpdmUgdnMuIEV4cG9uZW50aWFsCiAgKiBBUk1BIC8gQVJJTUEgbW9kZWxzCiAgKiBGaWx0ZXJzOiBLYWxtYW4sIEhvbHQtV2ludGVycwogICogU3RvY2hhc3RpYyBQcm9jZXNzZXMsIFJhbmRvbSBXYWxrLCBXaWVuZXIgUHJvY2VzcywgQnJvd25pYW4gTW90aW9uCiAgCgojIyMjIFdoYXQgaXMgYSBzdGF0aW9uYXJ5IHRpbWUgc2VyaWVzPwoKIkEgdGltZSBzZXJpZXMgaXMgc3RhdGlvbmFyeSBpZiBhIHNoaWZ0IGluIHRpbWUgZG9lc24ndCBjYXVzZSBhIGNoYW5nZSBpbiB0aGUgc2hhcGUgb2YgdGhlIGRpc3RyaWJ1dGlvbi4gQmFzaWMgcHJvcGVydGllcyBvZiB0aGUgZGlzdHJpYnV0aW9uIGxpa2UgdGhlIG1lYW4sIHZhcmlhbmNlLCBhbmQgY292YXJpYW5jZSBhcmUgY29uc3RhbnQgb3ZlciB0aW1lLiIKCiFbU3RhdGlvbmFyeSB2cyBTZWFzb25hbCBUaW1lIHNlcmllcyAtIGh0dHA6Ly93d3cuc3RhdGlzdGljc2hvd3RvLmNvbS9zdGF0aW9uYXJpdHkvXShpbWFnZTEucG5nKQoKCgojIyMjIERpZmZlcmVudCB0eXBlcyBvZiBzdGF0aW9uYXJpdHk6ICAKCiAgKiBTdHJpY3QgU3RhdGlvbmFyaXR5IC0gYmFzaWNhbGx5IGRvZXNuJ3QgZXhpc3QuLi4gZXNzZW50aWFsbHkgbm90aGluZyBjaGFuZ2VzIGJhc2VkIG9mZiBvZiB0aW1lLi4uCiAgKiBGaXJzdC1PcmRlciBTdGF0aW9uYXJpdHkgLSBtZWFucyBuZXZlciBjaGFuZ2UsIGJ1dCBvdGhlciBzdGF0aXN0aWNzIGNhbiBjaGFuZ2UKICAqIFNlY29uZC1PcmRlciAoV2VhaykgU3RhdGlvbmFyaXR5IC0gY29uc3RhbnQgbWVhbiwgdmFyaWFuY2UsIGFuZCBhdXRvY292YXJpYW5jZSB0aGF0IGRvZXNuJ3QgY2hhbmdlIHdpdGggdGltZSAKICAqIFRyZW5kLVN0YXRpb25hcnkgTW9kZWxzIC0gZmx1Y3R1YXRlIGFyb3VuZCBkZXRlcm1pbmlzdGljIHRyZW5kIChzZXJpZXMgbWVhbikuIFRyZW5kIGlzIHR5cGljYWxseSBsaW5lYXIgb3IgcXVhZHJhdGljLCBhbXBsaXR1ZGUgb2YgZmx1Y3R1YXRpb25zIHN0YXlzIGNvbnN0YW50CiAgKiBEaWZmZXJlbmNlLVN0YXRpb25hcnkgTW9kZWxzIC0gb25lIG9yIG1vcmUgZGlmZmVyZW5jZXMgb2YgdGhlIGRhdGEgbmVlZCB0byBiZSB0YWtlbiB0byBtYWtlIGl0IHN0YXRpb25hcnkKCgojIyMjIFdoeSBkbyB3ZSBjYXJlIGFib3V0IHN0YXRpb25hcml0eT8gIAoKTW9zdCBzdGF0aXN0aWNhbCBtb2RlbHMgdGhhdCBoYXZlIGJlZW4gZGV2ZWxvcGVkIGFzc3VtZSBzdGF0aW9uYXJpdHkuCgoKICAKYGBge3IsIHdhcm5pbmc9RkFMU0V9CnNldC5zZWVkKDEwMCkKCiMgbW9udGhseSBkYXRhLCBub3JtYWxseSBkaXN0cmlidXRlZAojIDEwMCBzdG9sdHpjb2luIGF2ZXJhZ2UgdHJhZGUgcHJpY2UgcGVyIG1vbnRoCiMgU3RhbmRhcmQgZGV2aWF0aW9uIG9mIDMKZGF0ZXMgPSBzZXEoZnJvbSA9IGx1YnJpZGF0ZTo6eW1kKCcyMDExLTA3LTAxJyksIAogICAgICAgICAgICB0byA9IGx1YnJpZGF0ZTo6eW1kKCcyMDE4LTA2LTMwJyksIAogICAgICAgICAgICBieSA9ICdtb250aHMnKQoKbiA9IGxlbmd0aChkYXRlcykKbXUgPSAxMDAKc2QgPSAzCgp0aW1lX3NlcmllcyA9IHRpYmJsZShkYXRlID0gZGF0ZXMsIHN0b2x0emNvaW4gPSBybm9ybShuID0gbiwgbWVhbiA9IG11LCBzZCA9IHNkKSwgCiAgICAgICAgICAgICAgICAgICAgIHR5cGUgPSAnc3RhdGlvbmFyeScsIHRyZW5kID0gJ25vbmUnKQoKZ2dwbG90KHRpbWVfc2VyaWVzLCBhZXMoeCA9IGRhdGUsIHkgPSBzdG9sdHpjb2luKSkgKyAKICBnZW9tX2xpbmUoKSArIGdlb21fc21vb3RoKG1ldGhvZCA9IGxtLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZvcm11bGEgPSB5IH4geCwgc2UgPSBGQUxTRSkgKyAKICBnZW9tX2hsaW5lKGFlcyh5aW50ZXJjZXB0ID0gbWVhbihzdG9sdHpjb2luKSksIGNvbG9yPSJyZWQiKQpgYGAKCgoKIyMjIyBQZXJmb3JtIHJlZ3Jlc3Npb24gYW5kIGxvb2sgYXQgdGhlIHJlc2lkdWFscwoKYGBge3J9CnRzX2ZpdCA9IGxtKHN0b2x0emNvaW4gfiBkYXRlLCBkYXRhID0gdGltZV9zZXJpZXMpCnRzX2ZpdApgYGAKCgpgYGB7cn0KcmVzID0gdGliYmxlKHJlc2lkdWFscyA9IHJlc2lkKHRzX2ZpdCkpCmdncGxvdChyZXMsIGFlcyh4ID0gMTpucm93KHJlcyksIHkgPSByZXNpZHVhbHMpKSArIGdlb21fcG9pbnQoKSArIGdlb21fcnVnKHNpZGVzID0gJ2wnKQpgYGAKCmBgYHtyfQpnZ3Bsb3QocmVzLCBhZXMoeCA9IHJlc2lkdWFscykpICsgZ2VvbV9oaXN0b2dyYW0oYWVzKHkgPSAuLmRlbnNpdHkuLiksIGJpbnMgPSA0MCkgKyBnZW9tX2RlbnNpdHkoY29sID0gJ3JlZCcsIHNpemUgPSAxLjUpCmBgYAoKCiMjIEFsd2F5cyB0ZXN0IGZvciBzdGF0aW9uYXJpdHkgZXZlbiBpZiBpdCBsb29rcyBzdGF0aW9uYXJ5CiMjIyMjIEF1Z21lbnRlZCBEaWNrZXktRnVsbGVyIGlzIGEgdHlwaWNhbCB0ZXN0IC0gdHNlcmllczo6YWRmLnRlc3QoKQogICogVGVzdHMgbnVsbCBoeXBvdGhlc2lzIHRoYXQgdGhlIGF1dG9yZWdyZXNzaXZlIG1vZGVsIGhhcyBhIHVuaXQgcm9vdAogICogU3RhdGlvbmFyaXR5IG1lYW5zIHJlamVjdGluZyB0aGUgbnVsbCBoeXBvdGhlc2lzCiAgKiBXZSB3aWxsIGxvb2sgZm9yIHAtdmFsdWUgPCAwLjA1CgoKYGBge3J9CiMgYWRmLnRlc3Qgb25seSBhY2NlcHRzIHVuaXZhcmlhdGUgdGltZSBzZXJpZXMgZGF0YSwgem9vIG9iamVjdCB3aWxsIGJlIHVzZWQKYWRmLnRlc3Qoem9vOjp6b28odGltZV9zZXJpZXMkc3RvbHR6Y29pbiwgdGltZV9zZXJpZXMkZGF0ZSkpCgpgYGAKCgoKIyMjIEFkZCBhIGxpbmVhciB0cmVuZCBmb3IgdGhlIGV4YW1wbGUKIyMjIyBTdG9sdHpjb2luIGhhcyBiZWVuIGdvaW5nIHVwIGJ5IDIgZWFjaCBtb250aAojIyMjIEtub3duIGFzICJ0cmVuZCBzdGF0aW9uYXJ5IgoKYGBge3J9CmdvaW5nX3VwX2J5ID0gMgp0cmVuZCA9IHNlcShmcm9tID0gMCwgCiAgICAgICAgICAgIHRvID0gKG4tMSkgKiBnb2luZ191cF9ieSwKICAgICAgICAgICAgYnkgPSBnb2luZ191cF9ieSkKCnRzX2xpbl90cmVuZCA9IHRpYmJsZShkYXRlID0gZGF0ZXMsIHN0b2x0emNvaW4gPSB0aW1lX3NlcmllcyRzdG9sdHpjb2luICsgdHJlbmQsIAogICAgICAgICAgICAgICAgICAgICB0eXBlID0gJ25vbi1zdGF0aW9uYXJ5JywgdHJlbmQgPSAnbGluZWFyJykKCmdncGxvdCh0c19saW5fdHJlbmQsIGFlcyh4ID0gZGF0ZSwgeSA9IHN0b2x0emNvaW4pKSArIAogIGdlb21fbGluZSgpICsgZ2VvbV9zbW9vdGgobWV0aG9kID0gbG0sIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgZm9ybXVsYSA9IHkgfiB4LCBzZSA9IEZBTFNFKSArIAogIGdlb21faGxpbmUoYWVzKHlpbnRlcmNlcHQgPSBtZWFuKHN0b2x0emNvaW4pKSwgY29sb3I9InJlZCIpCmBgYAoKCgojIyMjIFBlcmZvcm0gcmVncmVzc2lvbiBhbmQgbG9vayBhdCB0aGUgcmVzaWR1YWxzCgpgYGB7cn0KdHNfZml0ID0gbG0oc3RvbHR6Y29pbiB+IGRhdGUsIGRhdGEgPSB0c19saW5fdHJlbmQpCnRzX2ZpdApgYGAKCgpgYGB7cn0KcmVzID0gdGliYmxlKHJlc2lkdWFscyA9IHJlc2lkKHRzX2ZpdCkpCmdncGxvdChyZXMsIGFlcyh4ID0gMTpucm93KHJlcyksIHkgPSByZXNpZHVhbHMpKSArIGdlb21fcG9pbnQoKSArIGdlb21fcnVnKHNpZGVzID0gJ2wnKQpgYGAKCmBgYHtyfQpnZ3Bsb3QocmVzLCBhZXMoeCA9IHJlc2lkdWFscykpICsgZ2VvbV9oaXN0b2dyYW0oYWVzKHkgPSAuLmRlbnNpdHkuLiksIGJpbnMgPSA0MCkgKyBnZW9tX2RlbnNpdHkoY29sID0gJ3JlZCcsIHNpemUgPSAxLjUpCmBgYAoKCmBgYHtyfQphZGYudGVzdCh6b286Onpvbyh0c19saW5fdHJlbmQkc3RvbHR6Y29pbiwgdHNfbGluX3RyZW5kJGRhdGUpKQpgYGAKCgoKIyMjIEFkZCBhIHNlYXNvbmFsIGNvbXBvbmVudCB0byB0aGUgbGluZWFyIHRyZW5kCiMjIyMgU3RvbHR6Y29pbiBmbHVjdHVhdGVzIHRocm91Z2hvdXQgdGhlIHllYXIKCmBgYHtyfQpnb2luZ191cF9ieSA9IDIKdHJlbmQgPSBzZXEoZnJvbSA9IDAsIAogICAgICAgICAgICB0byA9IChuLTEpICogZ29pbmdfdXBfYnksCiAgICAgICAgICAgIGJ5ID0gZ29pbmdfdXBfYnkpCnNlYXNvbmFsID0gc2luKDE6bikgKiAxMAoKdHNfbGluX3RyZW5kX3NlYXNvbmFsID0gdGliYmxlKGRhdGUgPSBkYXRlcywgc3RvbHR6Y29pbiA9IHRpbWVfc2VyaWVzJHN0b2x0emNvaW4gKyB0cmVuZCArIHNlYXNvbmFsLCAKICAgICAgICAgICAgICAgICAgICAgdHlwZSA9ICdub24tc3RhdGlvbmFyeScsIHRyZW5kID0gJ2xpbmVhcicpCgpnZ3Bsb3QodHNfbGluX3RyZW5kX3NlYXNvbmFsLCBhZXMoeCA9IGRhdGUsIHkgPSBzdG9sdHpjb2luKSkgKyAKICBnZW9tX2xpbmUoKSArIGdlb21fc21vb3RoKG1ldGhvZCA9IGxtLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZvcm11bGEgPSB5IH4geCwgc2UgPSBGQUxTRSkgKyAKICBnZW9tX2hsaW5lKGFlcyh5aW50ZXJjZXB0ID0gbWVhbihzdG9sdHpjb2luKSksIGNvbG9yPSJyZWQiKQpgYGAKCgoKCiMjIyMgUGVyZm9ybSByZWdyZXNzaW9uIGFuZCBsb29rIGF0IHRoZSByZXNpZHVhbHMKCmBgYHtyfQp0c19maXQgPSBsbShzdG9sdHpjb2luIH4gZGF0ZSwgZGF0YSA9IHRzX2xpbl90cmVuZF9zZWFzb25hbCkKdHNfZml0CmBgYAoKCmBgYHtyfQpyZXMgPSB0aWJibGUocmVzaWR1YWxzID0gcmVzaWQodHNfZml0KSkKZ2dwbG90KHJlcywgYWVzKHggPSAxOm5yb3cocmVzKSwgeSA9IHJlc2lkdWFscykpICsgZ2VvbV9wb2ludCgpICsgZ2VvbV9ydWcoc2lkZXMgPSAnbCcpCmBgYAoKYGBge3J9CmdncGxvdChyZXMsIGFlcyh4ID0gcmVzaWR1YWxzKSkgKyBnZW9tX2hpc3RvZ3JhbShhZXMoeSA9IC4uZGVuc2l0eS4uKSwgYmlucyA9IDQwKSArIGdlb21fZGVuc2l0eShjb2wgPSAncmVkJywgc2l6ZSA9IDEuNSkKYGBgCgoKYGBge3J9CmFkZi50ZXN0KHpvbzo6em9vKHRzX2xpbl90cmVuZF9zZWFzb25hbCRzdG9sdHpjb2luLCB0c19saW5fdHJlbmRfc2Vhc29uYWwkZGF0ZSkpCmBgYAoKCgoKIyMjIEFkZCBhbiBleHBvbmVudGlhbCB0cmVuZCB0byB0aGUgc2Vhc29uYWwgY29tcG9uZW50CiMjIyMgU3RvbHR6Y29pbiBoYXMgaW5jcmVhc2VkIGluIHZvbGF0aWxpdHkgb3ZlciB0aGUgeWVhcnMKICAqIFN0b2x0emNvaW4gbm93IGhhczoKICAgICogQSBwb3NpdGl2ZSBhbmQgZ3Jvd2luZyBsaW5lYXIgdHJlbmQKICAgICogU2Vhc29uYWxpdHkKCmBgYHtyfQpnb2luZ191cF9ieSA9IDIKdHJlbmQgPSBzZXEoZnJvbSA9IDAsIAogICAgICAgICAgICB0byA9IChuLTEpICogZ29pbmdfdXBfYnksCiAgICAgICAgICAgIGJ5ID0gZ29pbmdfdXBfYnkpCnNlYXNvbmFsID0gc2luKDE6bikgKiAxOm4KCnRzX2xpbl90cmVuZF9zZWFzb25hbF9leHAgPSB0aWJibGUoZGF0ZSA9IGRhdGVzLCBzdG9sdHpjb2luID0gdGltZV9zZXJpZXMkc3RvbHR6Y29pbiArIHRyZW5kICsgc2Vhc29uYWwsIAogICAgICAgICAgICAgICAgICAgICB0eXBlID0gJ25vbi1zdGF0aW9uYXJ5JywgdHJlbmQgPSAnbGluZWFyJykKCmdncGxvdCh0c19saW5fdHJlbmRfc2Vhc29uYWxfZXhwLCBhZXMoeCA9IGRhdGUsIHkgPSBzdG9sdHpjb2luKSkgKyAKICBnZW9tX2xpbmUoKSArIGdlb21fc21vb3RoKG1ldGhvZCA9IGxtLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZvcm11bGEgPSB5IH4geCwgc2UgPSBGQUxTRSkgKyAKICBnZW9tX2hsaW5lKGFlcyh5aW50ZXJjZXB0ID0gbWVhbihzdG9sdHpjb2luKSksIGNvbG9yPSJyZWQiKQpgYGAKCgoKCiMjIyMgUGVyZm9ybSByZWdyZXNzaW9uIGFuZCBsb29rIGF0IHRoZSByZXNpZHVhbHMKCmBgYHtyfQp0c19maXQgPSBsbShzdG9sdHpjb2luIH4gZGF0ZSwgZGF0YSA9IHRzX2xpbl90cmVuZF9zZWFzb25hbF9leHApCnRzX2ZpdApgYGAKCgpgYGB7cn0KcmVzID0gdGliYmxlKHJlc2lkdWFscyA9IHJlc2lkKHRzX2ZpdCkpCmdncGxvdChyZXMsIGFlcyh4ID0gMTpucm93KHJlcyksIHkgPSByZXNpZHVhbHMpKSArIGdlb21fcG9pbnQoKSArIGdlb21fcnVnKHNpZGVzID0gJ2wnKQpgYGAKCmBgYHtyfQpnZ3Bsb3QocmVzLCBhZXMoeCA9IHJlc2lkdWFscykpICsgZ2VvbV9oaXN0b2dyYW0oYWVzKHkgPSAuLmRlbnNpdHkuLiksIGJpbnMgPSA0MCkgKyBnZW9tX2RlbnNpdHkoY29sID0gJ3JlZCcsIHNpemUgPSAxLjUpCmBgYAoKCmBgYHtyfQphZGYudGVzdCh6b286Onpvbyh0c19saW5fdHJlbmRfc2Vhc29uYWxfZXhwJHN0b2x0emNvaW4sIHRzX2xpbl90cmVuZF9zZWFzb25hbF9leHAkZGF0ZSkpCmBgYAoKCiMjIFIncyBkZWNvbXBvc2UgZnVuY3Rpb24gZG9lcyBhIGxvdCBmb3IgeW91CmBgYHtyfQpkMXkgPSAgYXMubnVtZXJpYyhsdWJyaWRhdGU6OnllYXIobWluKHRzX2xpbl90cmVuZF9zZWFzb25hbF9leHAkZGF0ZSkpKQpkMW0gPSBhcy5udW1lcmljKGx1YnJpZGF0ZTo6bW9udGgobWluKHRzX2xpbl90cmVuZF9zZWFzb25hbF9leHAkZGF0ZSkpKQpkMnkgPSAgYXMubnVtZXJpYyhsdWJyaWRhdGU6OnllYXIobWF4KHRzX2xpbl90cmVuZF9zZWFzb25hbF9leHAkZGF0ZSkpKQpkMm0gPSBhcy5udW1lcmljKGx1YnJpZGF0ZTo6bW9udGgobWF4KHRzX2xpbl90cmVuZF9zZWFzb25hbF9leHAkZGF0ZSkpKQoKeCA9IHRzKHRzX2xpbl90cmVuZF9zZWFzb25hbF9leHAkc3RvbHR6Y29pbiwgc3RhcnQ9YyhkMXksIGQxbSksIGVuZD1jKGQyeSwgZDJtKSwgZnJlcXVlbmN5PTEyKQoKcGxvdChkZWNvbXBvc2UoeCwgdHlwZSA9ICdhZGRpdGl2ZScpKQoKYGBgCgoKYGBge3J9CnBsb3QoZGVjb21wb3NlKHgsIHR5cGUgPSAnbXVsdGlwbGljYXRpdmUnKSkKYGBgCgoKIyMjIEV4YW1wbGUgZnJvbSBkZWNvbXBvc2UgZnVuY3Rpb24KIyMjIyByYW5kb20gcG9ydGlvbiBzaG91bGQganVzdCBsb29rIGxpa2Ugd2hpdGUgbm9pc2UgKG90aGVyd2lzZSB5b3UndmUgbGlrZWx5IG1pc3NlZCBzb21ldGhpbmcpCmBgYHtyfQpyZXF1aXJlKGdyYXBoaWNzKQoKbSA8LSBkZWNvbXBvc2UoY28yKQpwbG90KG0pCmBgYAoKCgpgYGB7ciwgd2FybmluZz1GQUxTRSwgbWVzc2FnZT1GQUxTRX0KZG93ID0gdHFfaW5kZXgoIkRPVyIpICU+JQogIHRxX2dldChnZXQgPSAic3RvY2sucHJpY2VzIikKaGVhZChkb3csIDEwKQpgYGAKCgpgYGB7ciB3YXJuaW5nPUZBTFNFLCBtZXNzYWdlPUZBTFNFfQpwID0gZ2dwbG90KGRvdywgYWVzKHggPSBkYXRlLCB5ID0gY2xvc2UsIGNvbCA9IHN5bWJvbCkpICsgCiAgZ2VvbV9saW5lKCkgKyAKICBnZ2hpZ2hsaWdodChtYXgoY2xvc2UpID4gMjUwKSArIAogIHRoZW1lX21pbmltYWwoKSArIAogIGZhY2V0X3dyYXAofiBzeW1ib2wpCnByaW50KHApCmBgYAoKCgojIyBJcyBCQSBhIHN0YXRpb25hcnkgdGltZSBzZXJpZXM/CgojIyMjIENhbiB5b3Ugc2ltcGx5IHpvb20gaW4gYW5kIGZpbmQgYSBwb3J0aW9uIHRoYXQgaXMgc3RhdGlvbmFyeT8KYGBge3J9CmJhID0gZG93ICU+JSAKICBmaWx0ZXIoc3ltYm9sID09ICdCQScpICU+JSAKICBzZWxlY3QoZGF0ZSwgY2xvc2UpCnBsb3RseTo6Z2dwbG90bHkoZ2dwbG90KGJhLCBhZXMoeCA9IGRhdGUsIHkgPSBjbG9zZSkpICsgZ2VvbV9saW5lKCkgKyBnZW9tX3Ntb290aChtZXRob2QgPSBsbSwgc2UgPSBGQUxTRSkpCmBgYAoKCmBgYHtyfQphZGYudGVzdCh6b286OnpvbyhiYSRjbG9zZSwgYmEkZGF0ZSkpCmBgYAoKCgpgYGB7cn0KdHNfZml0ID0gbG0oY2xvc2UgfiBkYXRlLCBkYXRhID0gYmEpCnRzX2ZpdApgYGAKCgpgYGB7cn0KcmVzID0gdGliYmxlKHJlc2lkdWFscyA9IHJlc2lkKHRzX2ZpdCkpCmdncGxvdChyZXMsIGFlcyh4ID0gMTpucm93KHJlcyksIHkgPSByZXNpZHVhbHMpKSArIGdlb21fcG9pbnQoKSArIGdlb21fcnVnKHNpZGVzID0gJ2wnKQpgYGAKCmBgYHtyfQpnZ3Bsb3QocmVzLCBhZXMoeCA9IHJlc2lkdWFscykpICsgZ2VvbV9oaXN0b2dyYW0oYWVzKHkgPSAuLmRlbnNpdHkuLiksIGJpbnMgPSA0MCkgKyBnZW9tX2RlbnNpdHkoY29sID0gJ3JlZCcsIHNpemUgPSAxLjUpCmBgYAoKCgoKCgojIyMjIEFyZSB0aGUgZGFpbHkgcmV0dXJucyBzdGF0aW9uYXJ5PwoKYGBge3J9CnJldCA9IGJhICU+JQogIHRxX3RyYW5zbXV0ZShjbG9zZSwgbXV0YXRlX2Z1biA9IGRhaWx5UmV0dXJuKQpnZ3Bsb3QocmV0LCBhZXMoeCA9IGRhdGUsIHkgPSBkYWlseS5yZXR1cm5zKSkgKyBnZW9tX2xpbmUoKSArIGdlb21fc21vb3RoKG1ldGhvZCA9IGxtLCBzZSA9IEZBTFNFKQpgYGAKCgpgYGB7cn0KYWRmLnRlc3Qoem9vOjp6b28ocmV0JGRhaWx5LnJldHVybnMsIHJldCRkYXRlKSkKYGBgCgoKYGBge3J9CnRzX2ZpdCA9IGxtKGRhaWx5LnJldHVybnMgfiBkYXRlLCBkYXRhID0gcmV0KQp0c19maXQKYGBgCgoKYGBge3J9CnJlcyA9IHRpYmJsZShyZXNpZHVhbHMgPSByZXNpZCh0c19maXQpKQpnZ3Bsb3QocmVzLCBhZXMoeCA9IDE6bnJvdyhyZXMpLCB5ID0gcmVzaWR1YWxzKSkgKyBnZW9tX3BvaW50KCkgKyBnZW9tX3J1ZyhzaWRlcyA9ICdsJykKYGBgCgpgYGB7cn0KZ2dwbG90KHJlcywgYWVzKHggPSByZXNpZHVhbHMpKSArIGdlb21faGlzdG9ncmFtKGFlcyh5ID0gLi5kZW5zaXR5Li4pLCBiaW5zID0gNDApICsgZ2VvbV9kZW5zaXR5KGNvbCA9ICdyZWQnLCBzaXplID0gMS41KQpgYGAKCgpgYGB7cn0KCmBgYAoK